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Abstract 

Background: A common survival strategy of microorganisms subjected to stress involves the generation 

of phenotypic heterogeneity in the isogenic microbial population enabling a subset of the population to survive 

under stress. In a recent study, a mycobacterial population of M. smegmatis was shown to develop phenotypic 

heterogeneity under nutrient depletion. The observed heterogeneity is in the form of a bimodal distribution of 

the expression levels of the Green Fluorescent Protein (GFP) as reporter with the gfp fused to the promoter of 

the rel gene. The stringent response pathway is initiated in the subpopulation with high rel activity. 

y—i , Results: In the present study, we characterise quantitatively the single cell promoter activity of the three 

key genes, namely, mprA, sigE and rel, in the stringent response pathway with gfp as the reporter. The origin of 

bimodality in the GFP distribution lies in two stable expression states, i.e., bistability. We develop a theoretical 

model to study the dynamics of the stringent response pathway. The model incorporates a recently proposed 

f*) ■ mechanism of bistability based on positive feedback and cell growth retardation due to protein synthesis. Based 

D ' on flow cytometry data, we establish that the distribution of GFP levels in the mycobacterial population at any 

Ph . point of time is a linear superposition of two invariant distributions, one Gaussian and the other lognormal, with 

~f^ ' only the coefficients in the linear combination depending on time. This allows us to use a binning algorithm 

and determine the time variation of the mean protein level, the fraction of cells in a subpopulation and also the 

i i ■ coefficient of variation, a measure of gene expression noise. 

Conclusions: The results of the theoretical model along with a comprehensive analysis of the flow cytometry 
^/^ , data provide definitive evidence for the coexistence of two subpopulations with overlapping protein distributions. 

.2" 

X> ! Background 

i__i , Microorganisms are subjected to a number of stresses during their lifetime. Examples of such stresses are: depletion 

of nutrients, environmental fluctuations, lack of oxygen, application of antibiotic drugs etc. Microorganisms take 

%jj recourse to a number of strategics for survival under stress and adapting to changed circumstances [HI]- A 

£-*. • prominent feature of such strategies is the generation of phenotypic heterogeneity in an isogenic microbial population. 

r- — , | The heterogeneity is advantageous as it gives rise to variant subpopulations which are better suited to persist under 

f^ ■ stress. Bistability refers to the appearance of two subpopulations with distinct phenotypic characteristics [5j[6]. In 

one of the subpopulations, the expression of appropriate stress response genes is initiated resulting in adaptation . 

£~^. ■ There are broadly two mechanisms for the generation of phenotypic heterogeneity [T1IB] ■ In "responsive switching" 

cells switch phenotypes in response to perturbations associated with stress. In the case of "spontaneous stochastic 

3 ■ switching" , transitions occur randomly between the phenotypes even in the absence of stress. Responsive switching 

may also have a stochastic component as fluctuations in the level of a key regulatory molecule can activate the 

switch once a threshold level is crossed [TJE1EM] • 

The pre-existing phenotypic heterogeneity, an example of the well-known "bet-hedging-strategy" , keeps the pop- 
ulation in readiness to deal with future calamities. Using a microfluidic device, Balaban ct al. [5] have demonstrated 
the existence of two distinct subpopulations, normal and pcrsister, in a growing colony of E. coli cells. The persistcr 
subpopulation constitutes a small fraction of the total cell population and is distinguished from the normal subpop- 
ulation by a reduced growth rate. Since killing by antibiotic drugs like ampicillin depends on the active growth of 
cell walls, the persister cells manage to survive when the total population is subjected to antibiotic treatment. The 
normal cells, with an enhanced growth rate are, however, unable to escape death. Once the antibiotic treatment 
is over, some surviving cells switch from the persister to the normal state so that normal population growth is 
resumed 0110] ■ A simple theoretical model involving transitions between the normal and persister phenotypes 
explains the major experimental observations well [9][lT]. In the case of environmental perturbations, Thattai and 
Oudenaarden [12] have shown through mathematical modelling that a dynamically heterogeneous bacterial pop- 
ulation can under certain circumstances achieve a higher net growth rate conferring a fitness advantage than a 
homogeneous one. Mathematical modelling further shows that responsive switching is favoured over spontaneous 
switching in the case of rapid environmental fluctuations whereas the reverse is true when environmental perturba- 
tions are infrequent [7] . Another theoretical prediction that cells may tune the switching rates between phenotypes 
to the frequency of environmental changes has been verified in an experiment by Acar et al. [13) involving an en- 
gineered strain of S. cerevisiae which can switch randomly between two phenotypes. The major feature of all such 
studies is the coexistence of two distinct subpopulations in an isogenic population and their interconversions in the 
presence/absence of stress. Bistability, i.e., the partitioning of a cell population into two distinct subpopulations 






has been experimentally observed in a number of cases [UGIM] . Some prominent examples include: lysis/lysogeny in 
bacteriophage A [14] , the activation of the lactose utilization pathway in E. coli [15| and the galactose utilization 
genetic circuit in S. cerevisiae |16j . competence development in B. subtilis [6l ll7[fT8j and the stringent response in 
mycobacteria [19]. 

The mycobacterial pathogen M. tuberculosis, the causative agent of tuberculosis, has remarkable resilience 
against various physiological and environmental stresses including that induced by drugs [20H22] . On tubercular 
infection, granulomas form in the host tissues enclosing the infected cells. Mycobacteria encounter a changed physical 
environment in the confined space of granulomas with a paucity of life-sustaining constituents like nutrients, oxygen 
and iron [231 HI]- The pathogens adapt to the stressed conditions and can survive over years in the so-called 
latent state. In vitro too, M. tuberculosis has been found to persist for years in the latent state characterised by 
the absence of active replication and metabolism |25j . Researchers have developed models simulating the possible 
environmental conditions in the granulomas. One such model is the adaptation to nutrient-depleted stationary phase 
[26] . The processes leading to the slowdown of replicative and metabolic activity constitute the stringent response. 
In mycobacteria, the expression of rel initiates the stringent response which leads to persistence. The importance of 
Rel arises from the fact that it synthesizes the stringent response regulator ppGpp (guanosine tetraphosphate) [27j 
and is essential for the long-term survival of M. tuberculosis under starvation |28] and for prolonged life of the bacilli 
in mice [25] . 

Key elements of the stringent response and the ability to survive over long periods of time under stress are shared 
between the mycobacterial species M. tuberculosis and M. smegmatis [30j . Recent experiments provide knowledge of 
the stress signaling pathway in mycobacteria linking polyphosphate (poly P), the two-component system MprAB, 
the alternate sigma factor SigE and Rel [31) . In an earlier study [19) . we investigated the dynamics of rel-gfp 
expression (gfp fused with rel promoter) in M. smegmatis grown upto the stationary phase with nutrient depletion 
serving as the source of stress. In a flow cytometry experiment, we obtained evidence of a bimodal distribution in 
GFP levels and suggested that positive feedback in the stringent response pathway and gene expression noise are 
responsible for the creation of phenotypic heterogeneity in the mycobacterial population in terms of the expression 
of rel-gfp. Positive feedback gives rise to bistability (5][6], i.e., two stable expression states corresponding to low 
and high GFP levels. We further demonstrated hysteresis, a feature of bistability, in rel-gfp expression. The 
mathematical model developed by us to study the dynamics of the stringent response pathway predicted bistability 
in a narrow parameter regime which, however, lacks experimental support. In general, to obtain bistability a gene 
circuit must have positive feedback and cooperativity in the regulation of gene expression. Recently, Tan et al. [32] 
have proposed a new mechanism by which bistability arises from a noncoopcrativc positive feedback circuit and 
circuit-induced growth retardation. The novel type of bistability was demonstrated in a synthetic gene circuit. 
The circuit, embedded in a host cell, consists of a single positive feedback loop in which the protein product X of 
a gene promotes its own synthesis in a noncooperative fashion. The protein decay rate has two components, the 
degradation rate and the dilution rate due to cell growth. In the circuit considered, production of X slows down 
cell growth so that at higher concentrations of X, the rate of dilution of X is reduced. This generates a second 
positive feedback loop since increased synthesis of X proteins results in faster accumulation of the proteins so that 
the protein concentration is higher. The combination of two positive feedback loops gives rise to bistability in the 
absence of cooperativity. A related study by Klumpp et al [33J has also suggested that cell growth inhibition by a 
protein results in positive feedback. 

Results 

In this paper, we develop a theoretical model incorporating the effect of growth retardation due to protein synthesis 
[321134] . We provide some preliminary experimental evidence in support of the possibility. In our earlier study [19] , 
bimodality in the rel-gfp expression levels was observed. As a control, GFP expression driven by the constitutive 
hsp60 promoter was monitored as a function of time. A single bright population was observed at different times of 
growth (Figure S4 of [19]). The unimodal rather than bimodal distribution ruled out the possibility that clumping 
of mycobacterial cells and cell-to-cell variation of plasmid copy number were responsible for the observed bimodal 
fluorescence intensity distribution of rel promoter driven GFP expression. In the present study, we perform flow 
cytometry experiments to monitor mprA-gfp and sigE-gfp expression levels. The distribution of GFP levels in each 
case is found to be bimodal. We determine the probability distributions of the two subpopulations associated with 
low and high expression levels at different time points in the three cases of mprA-gfp, sigE-gfp and rel-gfp expression. 
In each case, the total distribution is a linear combination of two invariant distributions with the coefficients in the 
linear combination depending on time. The results of hysteresis experiments are also reported. 
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Figure 1. Schematic diagram of the stringent response pathway in M. smegmatis activated under 
nutrient depletion. MprB-P and MprA-P are the phosphorylatcd forms of MprB and MprA respectively. Poly 
P serves as the phosphate donor in the conversion of MprB to MprB-P. 



Mathematical modeling of the stress response pathway 

Figure 1 shows a sketch of the important components of the stress response pathway in M. smegmatis subjected to 
nutrient depletion J191I31J . The operon mprAB consists of two genes mprA and mprB which encode the histidinc 
kinase sensor MprB and its partner the cytoplasmic response regulator MprA respectively. The protein pair responds 
to environmental stimuli by initiating adaptive transcriptional programs. Polyphosphate kinase 1 (PPK1) catalyses 
the synthesis of polyphosphate (poly P) which is a linear polymer composed of several orthophosphate residues. 
Mycobacteria possibly encounter a phosphate-limited environment in macrophages. Surcka ct al. [31] proposed 
that poly P could play a critical role under ATP depletion by providing phosphate for utilisation by MprAB. A 
recent experiment |34| on a population of M. tuberculosis has established that the MTB gene ppkl is significantly 
upregulated due to phosphate starvation resulting in the synthesis of inorganic poliphosphate (poly P). The two- 
component regulatory system SenX3-RegX3 is known to be activated on phosphate starvation in both M. smegmatis 
|35| and M. tuberculosis |34j . In the latter case, RegX3 has been shown to regulate the expression of ppkl, a 
feature expected to be shared by M. smegmatis. In both the mycobacterial populations, poly P regulates the 
stringent response via the mprA-sigE-rel pathway [31] . In our experiments, nutrient depletion possibly gives rise 
to phosphate starvation. On activation of the mprAB operon, MprB autophosphorylatcs itself with poly P serving 
as the phosphate donor [3~T]l3"4"]. The phosphorylatcd MprB-P phosphorylates MprA via phosphotransfcr reactions. 
There is also evidence that MprB functions as a MprA-P (phosphorylatcd MprA) phosphatase. MprA-P binds the 
promoter of the mprAB operon to initiate transcription. A positive feedback loop is functional in the signaling 
network as the production of MprA brings about further MprA synthesis. The mprAB operon has a basal level of 
gene expression independent of the operation of the positive feedback loop. Once the mprAB operon is activated, 
MprA-P regulates the transcription of the alternate sigma factor gene sigE, which in turn controls the transcription 
of rel. We construct a mathematical model to study the dynamics of the above signaling pathway. The new feature 
included in the model takes into account the possibility that the production of stress-induced proteins like MprA 
and MprB slows down cell growth. This effectively generates a positive feedback loop as explained in Refs. [521153] . 
Figure 2(a) shows the mean amount of GFP fluorescence in the total mycobacterial population as measured in a 
flow cytometry experiment [mprA promoter fused with gfp) versus time. Figure 2(a) shows the specific growth rate 
of the cell population versus time. The inset shows the experimental growth curve for the mycobacterial population. 
The growth was monitored by recording the absorbancc values at 600 nm spcctrophotomctrically (see Methods). 
The specific growth rate at time t is given by jvTtT - St wnere N(t) is the number of mycobacterial cells at time 
t. Nutrient depletion limits growth and proliferation and culminates in the activation of stress response genes. It 
appears that in many cases rapid growth and stress response are mutually exclusive so that the production of a 
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Figure 2. Growth retardation due to protein synthesis, (a) Mean amount of GFP fluorescence in the case 
of mprA promoter fused with gfp and (b) specific growth rate fj, of mycobacterial population versus time in hours 
(h). 



stress response protein gives rise to a slower growth rate [36] • The balance between the expressions of growth- 
related and stress-induced genes determines the cellular phenotype with respect to growth rate and stress response. 
Persister cells in both E. coli [H[10] and mycobacteria J21j[22] have slow growth rates. In the case of M. smegmatis, 
we have already established that the slower growing persister subpopulation has a higher level of Rel, the initiator 
of stringent response, as compared to the normal subpopulation [19j . The new addition to our mathematical 
model |19| involves nonlinear protein decay rates arising from cell growth retardation due to protein synthesis. We 
briefly discuss the possible origins of the nonlinearity and its mathematical form 32,33]. The temporal rate of 
change of protein concentration is a balance between two terms: rate of synthesis and rate of decay. The decay rate 
constant (j e ff) has two components: the dilution rate due to cell growth (/i) and the natural decay rate constant 
(7), i.e., Jeff — A* + 7 where fi is the specific growth rate. In many cases, the expression of a protein results in cell 
growth retardation J35K33] . The general form of the specific growth rate in such cases is given by 
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where x denotes the protein concentration and <f>,0 are appropriate parameters. In Ref. [32], the expression for 
(i (Eq. (1)) is arrived at in the following manner. The Monod model [37] takes into account the effect of resource or 
nutrient limitation on the growth of bacterial cell population. The rate of change in the number of bacterial cells is 
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where the specific growth rate jj, is given by 
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In (3), s is the nutrient concentration and k the half saturation constant for the specific nutrient. When s = k, 
the specific growth rate attains its half maximal value (/x TOax is the maximum value of specific growth rate). The 
metabolic burden of protein synthesis affecting the growth rate is modeled by reducing the nutrient amount s by e, 
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The magnitude of e is assumed to be small and proportional to the protein concentration x. Following the 
procedure outlined in the Supplementary Information of |32j . namely, applying Taylor's expansion to (4) and 
putting e = \x (A is a constant), one obtains the expression in Eq. (1) with = M "^ s and 6 = -^r . Thus, 
the decay rate of proteins has the form —j e ffX = — (7 + /z) x where \i is given by Eq. (1). There are alternative 
explanations for the origin of the nonlinear decay term, e.g., the synthesis of a protein may retard cell growth if it 
is toxic to the cell [33] . In the case of mycobacteria, there is some experimental evidence of cell growth retardation 
brought about by protein synthesis. The response regulator MprA has an essential role in the stringent response 
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Figure 3. Specific growth rate /i versus GFP fluorescence intensity xgfp fitted with an expression 
similar to that given in Eq. (1). The values of /4*£f and 9 GFP are fj, GFP = 0-94 and 6 GFP = 0.317. The 
data points correspond to the growth period of 16-23 hours. 



pathway leading to persistence of mycobacteria under nutrient deprivation. Inactivation of the regulator in an mprA 
insertion mutant resulted in reduced persistence in a murine model but the growth of the mutant was proved to be 
significantly higher than that observed in the cases of the wild- type species [351132] • 

Our experimental data (Figure 3) provide further support to the hypothesis that MprA synthesis leads to reduced 
specific growth rate. The data points represent GFP fluorescence intensity with gfp fused to the mprA promoter. 
The GFP acts as a reporter of the mprA promoter activity culminating in MprA (also MprB) synthesis. The data 
points shown in Figure 3 are those that correspond to the growth period of 16-23 hours in Figure 2. 

The data points are fitted by an expression similar to that in Eq. (1) with [i°^ a p = 0.94 and 8 GFP = 0.317. The 
differential equations describing the temporal rates of change of key protein concentrations in our model are described 
in the Additional File 1. Solving the equations, one finds the existence of bistability, i.e., two stable expression states 
in an extended parameter regime. Figures SI A-C (Additional File 1) show the plots for bistability and hysteresis 
for the proteins MprA, SigE and Rel versus the autophosphorylation rate. In the deterministic scenario and in the 
bistable regime, all the cells in a population are in the same steady state if exposed to the same environment and 
with the same initial state. The experimentally observed heterogeneity in a genetically identical cell population 
is a consequence of stochastic gene expression. The biochemical events involved in gene expression are inherently 
probabilistic [401141] in nature. The uncertainty introduces fluctuations (noise) around mean expression levels so 
that the single protein level of the deterministic case broadens into a distribution of levels. In the case of bistable 
gene expression, the distribution of protein levels in a population of cells is bimodal with two distinct peaks. 



Bimodal Expression of mprA, sigE and rel in M. smegmatis 

In the earlier study [19] , we investigated the dynamics of rel transcription in individual cells of M. smegmatis 
grown in nutrient medium up to the stationary phase, with nutrient depletion serving as the source of stress. We 
employed flow cytometry to monitor the dynamics of green fluorescent protein (GFP) expression in M. smegmatis 
harboring the rel promoter fused to gfp as a function of time. The experimental signature of bistaility lies in 
the coexistence of two subpopulations. We now extend the study to investigate the dynamics of mprA and sigE 
transcription in individual M. smegmatis cells in separate flow cytometry experiments. Figures 4(a) and (b) show 
the time course of mprA-GFP and sigE-GFP expressions respectively as monitored by flow cytometry. In both the 
cases, the distribution of GFP-expressing cells is bimodal indicating the existence of two distinct subpopulations. 
In each case, the cells initially belong to the subpopulation with low GFP expression. The fraction of cells with 
high GFP expression increases as a function of time. The two subpopulations with low and high GFP expression 
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Figure 4. Time course of (a) mprA-gfp and (b) sigE-gfp expression. M. smegmatis harboring the 
appropriate promoter construct was grown for different periods of time (indicated in hours (h)) and the specific 
promoter-driven expression of GFP was monitored by flow cytometry. With time, there is a gradual transition 
from the L to the H subpopulation. 

are designated as L and H subpopulations respectively. In the stationary phase, the majority of the cells belong to 
the H subpopulation. The presence of two distinct subpopulations confirms the theoretical prediction of bistability. 
We analysed the experimental data shown in Figure 4 and found that at any time point the distribution P{x, t) 
of GFP levels in a population of cells is a sum of two overlapping and time-independent distributions, one Gaussian 
(Pi (a;)) and the other lognormal (P 2 (x)), i.e., 



P(x,t)=C 1 (t)P 1 (x)+C 2 (t)P 2 (x) 



(5) 



The coefficients CVs (i=l, 2) depend on time whereas Pi(x) and P 2 (x) are time-independent. The general forms 
of P\{x) and P 2 (x) are, 



Px(x) 



P 2 (x) 



eM-(^) 2 ) 



woi 



exp{- 



1 l lnx — xp2 

2 V W 02 



) 2 ) 



X Wq 2 V 27T 



(6) 



(7) 



Figures S2(a) and (b) in Additional File 1 illustrate the typical forms of the Gaussian and lognormal dis- 
tributions. The Gaussian distribution has a symmetric form whereas the lognormal distribution is asymmetric 
and long-tailed. Figure 5 shows the experimental data for cell count versus GFP fluorescence intensity at se- 
lected time points in the cases when gfp is fused with mprA and sigE promoters in separate experiments. The 
dotted curves represent the individual terms in the r.h.s. of Eq. (5) and the solid curve denotes the linear com- 
bination P(x,t). The different parameters of P\{x) and P 2 (x) have the values xqi = 97.3366 (145.86181), Woi = 
103.0731 (154.67381), x 02 = 5.95526 (6.1171), w 02 = 0.17618(0.2509) when gfp is fused with mprA (sigE). The ratio 
of the coefficients, Ci(t)/C 2 (t), has the value listed by the side of each figure. Figure S3 displays a similar analysis 
of the experimental data when gfp is fused to the rel promoter. 

In the earlier study [19] , the total cell population was divided into L and H subpopulations depending on whether 
the measured GFP fluorescence intensity was less or greater than a threshold intensity. In the present study, we 
have obtained approximate analytic expressions for the distributions of GFP fluorescence intensity in the L and H 
subpopulations. The two distributions, Gaussian and lognormal, have overlaps in a range of fluorescence intensity 
values (Figure 5 and Figure S3 in Additional File 1). We next used the binning algorithm developed by Chang et 
al. |42| to partition the cells of the total population into two overlapping distributions, one Gaussian (Eq. (6)) and 
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Figure 5. Fitting of data with two distributions. Experimental data for cell count versus GFP fluorescence 
intensity at selected time points when gfp is fused with mprA and sigE promoters respectively. The solid curve 
represents P(x, t) in equation (5) and the dotted curves are the individual terms on the r.h.s. 



the other lognormal (Eq. (7)). At time t, let N(t) be the total number of cells. For each cell, the data Xj for the 
fluorescence intensity is used to calculate the ratios, 

Pi (X j ) + P 2 (Xj ) Pl[Xj) + Pi (Xj ) 

where P\{x) and Pz(x) are the distributions in Eqs. (6) and (7). A random number r is generated and the cell 
j is assigned to the L subpopulation if < r < gi(xj) , the cell belongs to the H subpopulation otherwise. Once 
the total population is partitioned into the L and H subpopulations, one can calculate the following quantities: 

«.W - *$ «-l,2) 

°?(<) = E(^(')-w(*)) 2 /w(*)-i) 

The indices i — 1, 2 correspond to the L and H subpopulations, w»(t) is the fraction of cells in the zth subpopula- 
tion at time t, Hi(i) is the mean fluorescence intensity for the ith. subpopulation and of (£) the associated variance. 
Figure 6 shows the results of the data analysis. Figure 6(a) shows the plots of mean GFP fluorescence level for 
the L subpopulation (basal level) versus time in the three cases of gfp fused with the promoters of mprA, sigE 
and rel respectively. Figure 6(b) displays the data for the fractions of cells, ui2(t) , versus time in the three cases 
and Figure 6(c) shows the transition rate versus time along with the coefficients of variation CV (CV= standard 
deviation/mean) of the protein levels in the L subpopulation versus time. 

Figure S4 (Additional File 1) shows the plots of mean GFP fluorescence level for the total population versus time 
in the three cases of gfp fused with the promoters of mprA, sigE and rel respectively. As in the case of the basal 
level versus time data (Figure 6(a)), the plots are sigmoidal in nature. We solved the differential equations of the 
theoretical model described in Additional File 1 and obtained the concentrations of MprA, MprB, SigE, MprA-P, 
MprB-P and GFP versus time. Some of these plots are shown in Figure S5 (Additional File 1) and reproduce the 
sigmoidal nature of the experimental plots. We note that the sigmoidal nature of the curves is obtained only when 
the non- linear nature of the degradation rate is taken into account. 

As we have already discussed, the distribution of GFP levels in the mycobacterial cell population is a linear 
combination of two invariant distributions, one Gaussian and the other lognormal, with only the coefficients in the 
linear combination dependent on time. Friedman et al. [43j have developed an analytical framework of stochastic 
gene expression and shown that the steady state distribution of protein levels is given by the gamma distribution. 
The theory was later extended to include the cases of transcriptional autoregulation as well as noise propagation in a 
simple genetic network. While experimental support for gamma distribution has been obtained earlier [44] , a recent 
exhaustive study |45| of the E. coli proteome and transcriptome with single-molecule sensitivity in single cells has 
established that the distributions of almost all the protein levels out of the 1018 proteins investigated, are well fitted 
by the gamma distribution in the steady state. The gamma distribution was found to give a better fit than the 
lognormal distribution for proteins with low expression levels and almost similar fits for proteins with high expression 
levels. We analysed our GFP expression data to compare the fits using lognormal and gamma distributions. For all 
the three sets of data (gfp fused with the promoters of mprA, sigE and rel), the lognormal and gamma distribution 
give similar fits at the different time points. Figure S6 (Additional File 1) shows a comparison of the fits for the 
case of gfp-mprA. The lognormal appears to give a somewhat better fit than the gamma distribution, specially at 
the tail ends. 

Hysteresis in gfp expression 

Some bistable systems exhibit hysteresis, i.e., the response of the system is history-dependent. In the earlier study, 
experimental evidence of hysteresis was obtained with gfp fused to the promoter of rel. The experimental procedure 
followed for the observation of hysteresis is as follows. In PPK-KO, the ppkl knockout mutant, the ppkl gene was 
introduced under the control of the tet promoter. We grew PPK-KO carrying the tetracycline-inducible ppkl and 
rel-gfp plasmid in medium with increasing concentration of tetracycline (inducer) . For each inducer concentration, 
the distribution of cells expressing gfp was analysed by flow cytometry in the stationary phase (steady state) and 
the mean GFP level was measured. A similar set of experiments was carried out for decreasing concentrations 
of tetracycline. In the present study, hysteresis experiments in the manner described above were carried out in 
the two cases of gfp fused to mprA and sigE promoters respectively. Figure 7 shows the hysteresis data (mean 



o MprA 



150 - 








- 


140 - 










130 - 










120 - 










110- 










100- 










90- 















10 


20 30 40 

Time (h) 


50 



80 
70 
60 
50 
40 
30 
20 
10- 





10 20 30 40 50 

Time (h) 







n 


10 


20 30 


40 50 


60 




B- 


















1 ' 


1 


























/ V 






m 




























a 














r 








/ /'\ 


\ r^ 


-■ 
















'(/) 


2- 






/ / \ 






CO 








' 1 \ 






h- 


0- 




"V 


1 x 




- 





















132 

130 = 
O 

128 § 
:d 
CL 

126 O 

124 w 



122 > 

o 



10 20 30 40 

Time (h) 



SigE 



SigE 




10 20 30 40 

Time (h) 




10 20 30 40 50 60 

Time (h) 




10 20 30 40 50 60 

Time (h) 



220- 








210- 








200- 








190- 








180- 








170- 








160- 








150- 






J^ 


140- 













10 


20 30 

Time (h) 



(a) 



Rel 



o Rel 



10 20 30 40 50 60 



40 50 60 




10 20 30 40 50 60 



Time (h) 

(b) 



IT 

5 2- 




10 20 30 40 50 60 

Time (h) 



(c) 



116 > 

o 



Figure 6. Analysis of the time course of gfp expression, (a) Mean protein level in L subpopulation (basal 
level) versus time in hours in the three cases of gfp fused with mprA, sigE and rel promoters respectively, (b) 
Fraction of cells W2(t) in the H subpopulation versus time in hours in the three cases, (c) Transition rate from the 
L to the H subpopulation and the CV of the protein levels in the L subpopulation versus time in hours in the 
three cases. The experimental data are analysed using binning algorithm to obtain the plots (a), (b) and (c). 
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Figure 7. Hysteresis in gfp expression. The gene gfp is fused with (a) mprA and (b) sigE promoter. Filled 
triangles and squares represent the experimental data of mean GFP fluorescence with increasing and decreasing 
concentrations of tetracycline inducer respectively. 

GFP fluorescence versus inducer concentration) in the two cases for increasing (branch going up) and decreasing 
(branch going down) inducer concentrations. The existence of two distinct branches is a confirmation of hysteresis 
in agreement with theoretical predictions (Figures SI A-C). Figure 8 shows the GFP distributions in the stationary 
phase for two sets of experiments with different histories, one in which the inducer concentration is increased from 
low to a specific value (indicated as "Low" in black) and the other in which the same inducer concentration is 
reached by decreasing the inducer concentration from a high value (indicated as "High" in red). The distributions 
show that two regions of monostability are separated by a region of bistability. In the cases of monostability, the 
distributions with different histories more or less coincide. In the region of bistability, the distributions are distinct 
indicating a persistent memory of initial conditions. 

Discussion 

The development of persistence in microbial populations subjected to stress has been investigated extensively in 
microorganisms like E. coli and mycobacteria [SlflO. 21, 22, 28. 2"5]. In an earlier study [15], we demonstrated the 
roles of positive feedback and gene expression noise in generating phenotypic heterogeneity in a population of M. 
smegmatis subjected to nutrient depletion. The heterogeneity was in terms of two distinct subpopulations designated 
as L and H subpopulations. The subpopulations corresponded to persister and non-persister cell populations with 
the stringent response being initiated in the former. In the present study, we have undertaken a comprehensive 
single cell analysis of the expression activity of the three key molecular players in the stringent response pathway, 
namely, MprA, SigE and Rel. This has been done by fusing gfp to the respective genes in separate experiments 
and monitoring the GFP levels in a population of cells via flow cytometry. The distribution has been found to be 
bimodal in each case. 

In our earlier study [T5], with only the positive autoregulation of the mprAB operon taken into account, bista- 
bility was obtained in a parameter regime with restricted experimental relevance. The inclusion of the effective 
positive feedback loop due to growth retardation by protein synthesis gives rise to a considerably more extended 
region of bistability in parameter space. The persister cells with high stringent response regulator levels are known 
to have slow growth rates [21], 22, 28 .29]. This is consistent with the view that stress response diverts resources from 
growth to stress- related functions resulting in the slow growth of stress-resistant cells [36] . Figures 2 and 3 provide 
experimental evidence that the mean intensity of GFP fluorescence monitoring mprA- gfp expression increases with 
time while the specific growth rate fx of the M. smegmatis population decreases in the same time interval. The 
reciprocal relationship between the two quantities is represented by an expression similar to that in Eq. (1). Since 
our knowledge of the detailed genetic circuitry involved in the stringent response is limited, we have not attempted 
to develop a model to explain the origin of cell growth retardation due to protein synthesis. Further experiments 
(e.g., sorting of the mycobacterial cell population into two subpopulations) are needed to provide conclusive evi- 
dence that increased protein synthesis retards cell growth. The stringent response pathway involving MprA and 
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Figure 8. Hysteresis via GFP distributions. The distributions in the stationary phase with two different 
histories (see text) when gfp is fused with (a) mprA and (b) sigE promoter. The specific inducer concentrations 
are mentioned with each plot. 



MprB is initiated when the mycobacterial population is subjected to stresses like nutrient depletion. There is 
now experimental evidence of complex transcriptional, translational, and posttranslational regulation of SigE in 
mycobacteria [46l449j . A double positive feedback loop arises due to the activation of transcription initiation of 
sigE by MprA-P and the activation of the transcription of the mprAB operon by the SigE-RNAP complex. Post- 
translational regulation of SigE is mediated by RseA, an anti-sigma factor. Barik et al. [49. have identified a novel 
positive feedback involving SigE and RseA which becomes functional under surface stress. More experiments need 
to be carried out to obtain insight on the intricate control mechanisms at work when mycobacteria are subjected 
to stresses like nutrient deprivation. This will lead to a better understanding of the major contributory factors 
towards the generation of phenotypic heterogeneity in mycobacterial populations subjected to stress. 



Conclusions 

In the present study, we have characterised quantatively the single cell promoter activity of three key genes in the 
stringent response pathway of the mycobacterial population M. smegmatis. Under nutrient depletion, a "responsive 
switching" occurs from the L to the H subpopulation with low and high expression levels respectively. A compre- 
hensive analysis of the flow cytometry data demonstrates the coexistence of two subpopulations with overlapping 
protein distributions. We have further established that the GFP distribution at any time point is a linear super- 
position of a Gaussian and a lognormal distribution. The coefficients in the linear combination depend on time 
whereas the component distributions arc time-invariant. The Gaussian and lognormal distributions describe the 
distribution of protein levels in the L and H subpopulations respectively. The two distributions overlap in a range 
of GFP fluorescence intensity values. We also find that the experimental data for the H subpopulation can be fitted 
very well by the gamma distribution though the lognormal distribution gives a slightly better fit. In the case of 
skewed positive data sets, the two distributions are often interchangeable |50j . An analytical framework similar 
to that in Ref. [13] is, however, yet to be developed for the mycobacterial stringent response pathway studied in 
the paper. The major components in the pathway are the two-component system mprAB and multiple positive 
feedback loops. The two-component system is known to promote robust input-output relations |51j and persistence 
of gene expression states [52] which may partly explain the good fitting of the experimental data by well-known 
distributions. Further quantitative measurements combined with appropriate stochastic modeling are needed to 
characterise the experimentally observed subpopulations more uniquely. We used the binning algorithm developed 
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in [42] to partition the experimental cell population into the L and H subpopulations. This enabled us to compute 
quantities like the mean protein level in the L subpopulation, the fraction of cells in the H subpopulation and the 
CV of GFP levels in the L subpopulation as a function of time. The picture that emerges from the analysis of 
experimental data is that of bistability, i.e., the coexistence of two distinct subpopulations and stochastic transitions 
between the subpopulations resulting in the time evolution of the fraction of cells in the H subpopulation. As pointed 
out in the earlier study [19] , the rate of transition to the H subpopulation and the CV of the L subpopulation levels 
attain their maximum values around the same time point (Figure 5(c)) indicating the role of gene expression noise 
in bringing about the transition from the L to the H subpopulation. We have not attempted to develop theoretical 
models describing the time evolution of the relative weights, tu^s (i = 1,2), of the two subpopulations (Eq. (9)). 
A simple model of two interacting and evolving subpopulations with linear first order kinetics |12j , cannot explain 
the sigmoidal nature of the time evolution. A model with nonlinear growth kinetics has been proposed in [42) but 
lacking definitive knowledge on the origin of nonlinearity in the growth of mycobacterial subpopulations we defer 
the task of model building to a future publication. 

Methods 

Strains 

M. smegmatis mc 2 155 was grown routinely in Middle Brook (MB) 7H9 broth (BD Biosciences) medium supple- 
mented with 2% glucose and 0.05% Twcen 80. 

Construction of plasmids for fluorescence measurements 

The mprAB promoter was amplified from the genomic DNA of M. smegmatis using the sense and antisense primers, 
5'-AAGGTACCGCGCAACACCACAAAAAGCG-3' and 5'-TAGGATCCAGTTTTGACTCACTATCTGAG-3' 
respectively and cloned into the promoter-less replicative gfp vector pFPV27 between the Kpnl and BamHI sites 
(in bold). The sigE and rel promoters fused to gfp have been described earlier [191131] . The resulting plasmids 
were electroporatcd into M. smegmatis mc 2 155 for further study. For the study of hysteresis, expression of ppkl 
under a tetracycline-inducible promoter in an M. smegmatis strain inactivated in the ppkl gene (PPK-KO), has 
been described earlier [19] . 

FACS analysis 

M. smegmatis cells expressing different promoters fused to GFP were grown in medium supplemented with kanamycin 
(25 fig /ml) and analysed at different points of time on a FACS Caliber (BD Biosciences) flow cytometer as described 
earlier [19"] . Briefly, cells were washed, resuspended in PBS and fluorescence intensity of 20,000 events was measured. 
The data was analyzed using Cell Quest Pro (BD Biosciences) and WINMIDI software. The flow cytometry data 
is represented in histogram plots where the x-axis is a measure of fluorescence intensity and the y-axis represents 
the number of events. 

Measurement of growth rate 

M. smegmatis expressing promoter- g/p fusion constructs were grown in Middle Brook (MB) 7H9 broth supplemented 
with glucose and Twecn 80, and kanamycin (25 fig /ml). Growth at different time points was measured by recording 
absorbance values at 600 nm (a value of 1 ODgoo is equal to 10 8 cells or 200 fig dry weight of cells). A growth curve 
was generated by plotting absorbance values against time (inset of Figure 2). The specific growth rate fi (Eq. (2)) 
at different time points is determined by taking derivatives of the growth curve at the different time points (Figure 
2). 

Authors contributions 

IB, JB and MK conceptualised, supervised and coordinated the study. KS, JB and MK carried out the experiments. 
SG, BG and IB developed the theoretical model, performed the data analysis and interpreted the data. IB drafted 
the manuscript. All authors read and approved of the final version. 



12 



Acknowledgements 

IB, JB and MK thank M. Thattai for some useful discussions. This work was supported in part by a grant from 
the Department of Biotechnology, Government of India to MK. SG is supported by CSIR, India, under Grant No. 
09/015(0361)/2009-EMR-I. 



13 



References 

1. Vccning JW, Smits WK, Kuipcrs OP: Bistability, epigenetics and bet-hedging in bacteria. Annu Rev 
Microbiol 2008, 62:193-210. 

2. Davidson CJ, Surette MG: Individuality in bacteria. Annu Rev Genet 2008, 42:11.1-11.6. 

3. Fraser D, Kacrn M: A chance at survival: gene expression noise and phenotypic diversification 
strategies. Mol Microbiol 2009, 71:1333-1340. 

4. Locke JCW, Elowitz MB: Using movies to analyse gene circuit dynamics in single cells. Nat Rev 

Microbiol 2009, 7:383392. 

5. Dubnau D, Losick R: Bistability in bacteria. Mol Microbiol 2006, 61:564-572. 

6. Smits WK, Kuipcrs OP, Vccning JW: Phenotypic variation in bacteria: the role of feedback regu- 
lation. Nat Rev Microbiol 2006, 4:259-271. 

7. Kussell E, Leibler S: Phenotypic diversity, population growth and information in fluctuating en- 
vironments. Science 2005, 309:2075-2078. 

8. Lcisncr M, Stingl K, Frey E, Maicr B: Stochastic switching to competence. Curr Opin Microbiol 2008, 
11:553-559. 

9. Balaban NQ, Merrin J, Chait R, Kowalik L, Leibler S: Bacterial persistence as a phenotypic switch. 

Science 2008, 305:1622-1625. 

10. Gefen O, Balaban NQ: The importance of being persistent: heterogeneity of bacterial populations 
under antibiotic stress. FEMS Microbiol Rev 2009, 33:704-717. 

11. Dhar N, McKinney JD: Microbial phenotypic heterogeneity and antibiotic tolerance. Curr Opin 
Microbiol 2007, 10:30-38. 

12. Thattai M, van Oudenaarden A: Stochastic gene expression in fluctuating environments. Genetics 
2004, 167:523-530. 

13. Acar M, Mettetal JT, van Oudenaarden A: Stochastic switching as a survival strategy in fluctuating 
environments. Nat Genet 2008, 40:471-479. 

14. Arkin A, Ross J, McAdams HH: Stochastic kinetic analysis of developmental pathway bifurcation 
in phage A-infected E. coli cells. Genetics 1998, 149:1633-1648. 

15. Ozbudak EM, Thattai M, Lim HN, Shraiman BI, van Oudenaarden A: Multistability in the lactose 
utilization network of E. coli. Nature 2004, 427:737-740. 

16. Acar M, Becskei M, van Oudenaarden A: Enhancement of cellular memory by reducing stochastic 
transitions. Nature 2004, 435:228-232. 

17. Sel G, Garcia-Ojalvo J, Liberman LM, Elowitz M: An excitable gene regulatory circuit induces tran- 
sient cellular differentiation. Nature 2006, 440:545-550. 

18. Maamar H, Raj A, Dubnau D: Noise in gene expression determines cell fate in B. subtilis. Science 
2007, 317:526-529. 

19. Sureka K, Ghosh B, Dasgupta A, Basu J, Kundu M, Bose I: Positive feedback and noise activate the 
stringent response regulator Rel in mycobacteria. PLoS ONE 2008, 3:c 1771. 

20. McCunc RM, Fcldman FM, Lambert HP, McDcrmont W: Microbial persistence. II. Characteristics of 
the sterile state of tubercle bacilli. J Exp Med 1966, 123:445-468. 

21. Lewis K: Persister cells, dormancy and infections disease. Nat Rev Microbiol 2007, 5:48-56. 

22. Young D, Stark J, Kirschner D: Systems biology of persistent infection: tuberculosis as a case 
study. Nat Rev Microbiol 2008, 6:520-528. 

14 



23. Bishai W: Lipid lunch for persistent pathogen. Nature 2000, 406:683685. 

24. De Voss JJ, Rutter K, Schrocdcr BG, Barry III CE: Iron acquisition and metabolism by mycobacteria. 

J Bacteriol 1999, 181:44434451. 

25. Nyka W: Studies on the effect of starvation on mycobacteria. Infect Immun 1974, 9:843850. 

26. Smeuldcrs MJ, Keer J, Speight RA, Williams HD: Adaptation of Mycobacterium smegmatis to sta- 
tionary phase. J Bacteriol 1999, 181:270-283. 

27. Braeken K, Moris M, Daniels R, Vanderleydcii J, Michiels J: New horizons for (p)ppGpp in bacterial 
and plant physiology. Trends Microbiol 2006, 14:45-54. 

28. Primm TP, Andersen SJ, Mizrahi V, Avarbock D, Rubin H, et al: The stringent response of Mycobac- 
terium tuberculosis is required for long-term survival. J Bacteriol 2000, 182:4889-4898. 

29. Dahl JL, Kraus CN, Boshoff HIM, Doan B, Foley K, ct al: The role of RelMtb-mediated adaptation 
to stationary phase in long-term persistence of Mycobacterium tuberculosis in mice. Proc Natl 
Acad Sci USA 2003, 100:10026-10031. 

30. Ojha AK, Mukhcrjee TK, Chatterji D: High intracellular level of guanosine tetraphosphate in My- 
cobacterium smegmatis changes the morphology of the bacterium. Infect Immun 2000, 68:4084 
4091. 

31. Surcka K, Dcy S, Datta P, Singh AK, Dasgupta A, et al: Polyphosphate kinase is involved in stress- 
induced mprAB-sigE-rel signaling in mycobacteria. Mol Microbiol 2007, 65:261276. 

32. Tan C, Marguet P, You L: Emergent bistability by a growth-modulating positive feedback circuit. 

Nat Chem Biol 2009, 5:842-848. 

33. Klumpp S, Zhang Z, Hwa T: Growth rate-dependent global effects on gene expression in bacteria. 
Cell 2009, 139:1366-1375. 

34. Glover RT, Kriakov J, Garforth SJ, Baughn AD, Jacobs WRJ: The two-component regulatory system 
senX3-regX3 regulates phosphate-dependent gene expression in Mycobacterium smegmatis. J 

Bacteriol 2007, 189:5495-5503. 

35. Rifat D, Bishai WR, Karakousis PC: Phosphate Depletion: A Novel Trigger for Mycobacterium 
tuberculosis Persistence. J Infect Dis 2009, 200:1126-1135. 

36. Lpcz-Maury L, Marguerat S, Bhler J: Tuning gene expression to changing environments: from rapid 
responses to evolutionary adaptation. Nat Rev Genet 2008, 9:583-593. 

37. Monod J: The growth of bacterial cultures. Annu Rev Microbiol 1949, 3:371-394. 

38. Zhart TC, Deretic V: Mycobacterium tuberculosis signal transduction system required for per- 
sistent infections. Proc natl Acad Sci USA 2001, 98:12706-12711. 

39. Pang X, Phong V, Thomas F B, Salccna G, et al: Evidence for complex interactions of stress- 
associated regulons in an mprAB deletion mutant of Mycobacterium tuberculosis. Microbiology 
2007, 153:1229-1242. 

40. Kaern M, Elston TC, Blake WJ, Collins JJ: Stochasticity in gene expression: from theories to phe- 

notypes. Nat Rev Genet 2005, 6:451464. 

41. Raj A, van Oudcnaarden A: Nature, Nurture, or Chance: Stochastic Gene Expression and Its 
Consequences. Cell 2008, 135:216-226. 

42. Chang HH, Hemberg M, Barahona M, Ingber DE, Huang S: Transcriptome-wide noise controls lineage 
choice in mammalian progenitor cells. Nature 2008, 453:544-547. 

43. Friedman N, Cai L, Xie XS: Linking Stochastic Dynamics to Population Distribution: An Analyt- 
ical Framework of Gene Expression. Phys. Rev. Lett. 2006, 97:168302(1)-168302(4). 

15 



44. Cai L, Friedman N, S XX: Stochastic protein expression in individual cells at the single molecule 

level. Nature 2006, 440:358-362. 

45. Taniguchi Y. Choi PJ, Li GW, Chen H, Babu M, et al: Quantifying E. coli Proteome and Transcrip- 
tome with Single-Molecule Sensitivity in Single Cells. Science 2010, 329 (5991):533-538. 

46. Don V, Rodrigue S, Dainese E, Pal G, Gaudreau L, Manganelli R, Provvedi R: Evidence of Complex Tran- 
scriptional, Translational, and Posttranslational Regulation of the Extracytoplasmic Function 
Sigma Factor a E in Mycobacterium tuberculosis. J Bacteriol 2008, 190:5963-5971. 

47. Rodrigue S, Provvedi R, Jacques PE, Gaudreau L, Manganelli R: The sigma factors of Mycobacterium 
tuberculosis. FEMS Microbiol Rev 2006, 30:926941. 

48. He H, Hovey R, Kane J, Singh V, Zahrt TC: MprAB is a stress-responsive two-component system that 
directly regulates expression of sigma factors SigB and SigE in Mycobacterium tuberculosis. J 

Bacteriol 2006, 188:2134-2143. 

49. Barik S, Sureka K, Mukherjee P, Basu J, Kundu M: RseA, the SigE specific anti-sigma factor of My- 
cobacterium tuberculosis, is inactivated by phosphorylation-dependent ClpClP2 proteolysis. 

Mol Microbiol 2009, 75:592606. 

50. Wiens BL: When log-normal and gamma models give different results: a case study. The American 
Statistician 1999, 53(2):89-93. 

51. Shinar G, R M, Martinez MR, U A: Input-output robustness in simple bacterial signaling systems. 
PNAS 2007, 104:19931-35. 

52. Kato A, Mitrophanov AY, Groisman EA: A connector of two-component regulatory systems pro- 
motes signal amplification and persistence of expression. PNAS 2007, 104:12063-68. 



16 



Supplementary Information 

Mathematial model 

The reaction scheme describing the processes shown in Figure 1 is given by 

B 2 +B 2 *± B 2 - B 2 A Bi + B l (10) 

k 2 

k 5 

h 

ks 

" : » 

kd 

G AB + A 1 -A A 2 +B 2 (14) 

In the equations, A\(A 2 ) represents the phosphorylated (unphosphorylatcd) form of MprA and B\(B 2 ) denotes 
the phosphorylated (unphosphorylatcd) form of MprB. The inactive and active states of the mprAB operon are 
represented by Gab and Gab* respectively. In the inactive state, MprA and MprB proteins are synthesized at 
a basal rate s and in the active state protein production occurs at an enhanced rate /3. Eq. (1) describes the 
autophosphorylation reaction of MprB with the poly P chain serving as a source of phosphate groups. Eq. (2) 
describes the transfer of the phosphate group from the phosphorylated MprB to MprA. Eq. (3) corresponds to 
dcphosphorylation of phosphorylated MprA by unphosphorylated MprB which thus acts as a phosphatase (in 
the earlier study pQ, phosphorylated MprB was assumed to act as a phosphatase which is not consistent with 
experimental evidence). Eqs. (4) and (5) describe activation of the mprAB operon by phosphorylated MprA and 
basal expression of the operon respectively. Refs. [5]-3] provide experimental justification for the reaction scheme 
shown in Eqs. (l)-(5). Using standard mass action kinetics, we write down the rate equations for the concentration 
of each of the key molecular species participating in the biochemical events. The equations are: 



A 2 +B x ^ At-Bi -^> A x +B 2 (11) 



A x + B 2 ^ A x - B 2 -^> A 2 + B 2 (12) 



Gab + A x ^ G AB -A A 2 + B 2 (13) 



|il = k 6 [A 2 B 1 ] k 7 [A x }[B 2 ] + fc 8 L4i - B 2 ] 7^] - jj^j 



(15) 



^T - • + ^ + ^-^-^*] + ^-^-^-IT& (16) 



d[B 



il = k 3 [B 2 - B 2 ] - fc 4 [A 2 pi] + k 5 [A 2 - B x ] - 7 [B X ] - — ffi 



dt ° l ' ZJ -i"-H— j ■ -r- —J n~v l + ^[Bi] 



(17) 



-k 2 [B 2 - B 2 ] + k x [B 2 } 2 - k 3 [B 2 - B 2 ] (19) 

*4[A a ][Bi] - h[A 2 - Bi] - k 6 [A 2 - Bi] (20) 

= k 7 [A x ][B 2 } - (fc 8 + k 9 )[A x - B 2 ] (21) 
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d[B 2 - B 2 


1 


dt 




d[A 2 - Bx] 




dt 




d[A x ~ 


B 2 ] 



d[SigE\ [A^/k' 

-*- =Sl + h i + ^j/y " Sl [SigE] (22) 

^[GFP] JStgEj/k^ 

—*T = S2+k 1 + [SigE\/k» &2 [GFP] (23) 

Eq. (13) represents SigE synthesis due to transcriptional activation of the sigE gene by phosphorylated MprA-P. 
Eq. (14) describes GFP production due to the activation of the rel promoter by SigE. The rate constants 7, 71, S\ 
and 8 2 are the degradation rate constants. The last terms in Eqs. (6)-(9) represent the nonlinear decay rates the 
genesis of which is explained in the main text (see Eq. 4) [5]. Eqs. (6)-(14) correspond to the case where gfp is fused 
to the rel promoter. In the other cases when gfp is fused to the mprA or sigE promoter, appropriate modifications 
in the set of equations are required. 

The steady state solution of Eqs. (6)-(14) is obtained by setting all the rates of change to be zero. In the case 
of bistability, there are three steady state solutions, two stable and one unstable [6H8]. In the steady state, one has 
to solve the following set of coupled nonlinear algebraic equations: 

aMlBi] a2 [Ai\[B2] l[A{\ 1 f^ M] = (24) 



where, 



a[B 2 } 2 - ai [A 2 ] [B{\ - 7 [Bi] - 1 f^ 1 = (26) 

- + ^TTjim - am ' + a ' lA ' ]m - 7[al - TT^rn = < 27 > 

°> + ^TTJ&EW- SAGFP, '° (29) 

hh k^ke k 7 k 9 k d 

a = u i u ' ai = I 1 i ' a2 = 1 — TIT' k = T~ ( 30 ) 

ki + k 3 k 5 + k e k$ + kg k a 

The solutions of Eqs. (15)-(20) are obtained with the help of Mathematica. Figures SI A-C show the steady 
state solutions generated by varying the parameter a (associated with the autophosphorylation of MprB). The 
parameters have values:^ = 2.4, a 2 = 2.8,7 = 0.1, s = 0.14, (3 = 4, k = 1,71 = l,<f> = 0.5, 61 = l,9 2 ~ 10, Si = 
0.02, j3 1 = 4, k' = 10, 61 = 1, s 2 = 0.12, /? 2 = 4, k" = 2,6 2 = 0.1 in appropriate units. 

In each of the Figures SI A-C, the solid and dotted branches represent stable and unstable steady states 
respectively. Bistability is obtained over a wide range of parameter values due to the inclusion of the non-linear 
decay term in Eqs. (6)-(9). In the hysteresis experiments, the inducer tetracycline was used to control the level 
of PPK1 and therefore the synthesis of the poly P chain. Since the latter acts as the source of phosphate groups 
for the autophosphorylation of MprB (Eq. (1)), the rate constant k, in Eq. (1) is effectively proportional to the 
inducer (or the PPK1) concentration. As the parameter a (Eq. (21)) includes the rate constant k\ , a varying 
inducer concentration is equivalent to varying the parameter a. There is some experimental evidence that MprA-P 
regulates the expression of the mprAB operon in the form of dimers [9|. Inclusion of this feature in our model 
makes the bistable behaviour more prominent. 
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Figure SI: Bistability and hysteresis in the deterministic model. Steady state concentrations of MprA, SigE and 
Rel versus the parameter a (Eq. (21)). 
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Figure S2: (a) Gaussian and (b) lognormal distributions which describe the distribution of GFP leads in the L 
and H subpopulations respectively. 
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Figure S3: Experimental data for cell count versus GFP fluorescence intensity at selected time points when gfp is 
fused with rel promoter. The solid curve represents P(x, t) in Eq. (5) and the dotted curves are the individual terms 
on the r.h.s. The different parameters of Pi (a;) and ^(^O have the values Xqi = 157.14748, w i = 150.43575, x 02 = 
6.10036, u>02 = 0.1847 when gfp is fused with rel. 
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Figure S4: Mean GFP fluorescence level for the total population versus time in the three cases of gfp fused with 
the promoters of (a) mprA, (b) sigE and (c) rel. 
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Figure S5: Concentrations of MprA, SigE and GFP (in arbitrary units) versus time. The values of the concentra- 
tions are obtained by solving Eqs. (6)-(14) in Text SI. 
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Figure S6: Comparison of fits of experimental data for cell count versus GFP fluorescence intensity at selected 
time points when gfp is fused with mprA promoter, with lognormal (Eq. (7)) and gamma distributions. The gamma 

distribution has the form P{x) = - — h ^r^ ~ : where the parameters a and b have the values a = 33.246, b = If. 59 
and r(a) is the gamma function. 
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